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We study dipole oscillations in a general fermionic mixture: starting from the Boltzmann equation, 
we classify the different solutions in the parameter space through the number of real eigenvalues of 
the small oscillations matrix. We discuss how this number can be computed using the Sturm algo- 
rithm and its relation with the properties of the Laplace transform of the experimental quantities. 
After considering two components in harmonic potentials having different trapping frequencies, we 
study dipole oscillations in three-component mixtures. Explicit computations are done for realis- 
tic experimental setups using the classical Boltzmann equation without intra-species interactions. 
A brief discussion of the application of this classification to general collective oscillations is also 
presented. 

PACS numbers: 67.85.Lm, 51.10.+y 



I. INTRODUCTION 

The study of collective modes is a major tool to charac- 
terize and unravel the effects of interparticle interactions 
in a broad range of physical systems and in particular 
in experiments of ultracold atoms, where the scattering 
length can be tuned through Feshbach resonances [l|, [|J ■ 
For ultracold bosons the collective oscillation frequen- 
cies in single-component gases depend on the interac- 
tion energy and in multi-component gases they depend 
on both intra- and inter-species interactions; for ultra- 
cold fcrmions the s-wave intra-species scattering length 
vanishes due to Pauli principle and the collective modes 
depend on the interactions between the different compo- 
nents, so that their study is crucial in multi-component 
fermionic mixtures. 

In the last decade, the impressive experimentalpro- 
gresses in trapping and controlling Fermi gases 0, 3 al- 
lowed their study in highly controllable setups: it became 
possible to lower the temperature to reveal degeneracy ef- 
fects [f|, tune the interactions to explore the BEC-BCS 
crossover superimpose optical lattices fl3rJl6| and 

polarize two-component Fermi mixtures [13, Il8| . Collec- 
tive modes in two-com pon ent mixtures have also been 
experimentally studied [19l - l24| . 

Among these progresses, it recently became ex- 
perimentally manageable to study three-component 
fermionic mixtures of 6 Li [H[ and of 6 Li - 40 K (Si^. 
In 2") a degenerate Fermi gas consisting of three different 
hyperfme states of 6 Li was created: the three scattering 
lengths eti2, ai3 and 023 are different and tunable; the 
collisional stability of the gas has been studied, showing 
rather long lifetimes when one scattering length is close 
to the unitary limit and the other two are small 25]. In 
[2r3 | the interspecies Feshbach resonances of a mixture of 
6 Li and 40 K were studied: the collisional stability has 



been recently determined 27jj , showing that a small sam- 
ple of 40 K immersed in a two-component 6 Li mixture is 



stable with very low loss rates for large negative scatter- 
ing lengths between the two 6 Li. 

These systems are particularly promising for several 
reasons: they offer a way to investigate new strongly 
interacting Fermi systems; thanks to their enlarged pa- 
rameter space, they could present novel phases related 
to the formation of unconventional pairing and stabil- 
ity of trions |28i -T2 : moreover, the possibility of having 
many components with tunable interactions can be used 
to study mechanisms of color superconductivity [33j and, 
in perspective, to simulate parts of the QCD phase dia- 
gram [34| . 

A natural way to characterize multi-component 
fermionic mixtures is provided by the study of their 
modes. For two-component Fermi mixtures, collective 
oscillations have been investigated in several situations 
of experimental relevance (23l. l35TJ42l | (see more references 
in the reviews and their behaviour studied when 

different parameters are varied: temperature, scattering 
length, polarization (i.e., the relative number of atoms of 
the components) and mass ratio. The variety of possi- 
ble dynamical regimes grows when the number of com- 
ponents increases: quantized vortices of three-component 
Fermi mixtures with attractive interactions have been re- 
cently studied (43j . 

The goal of the present paper is to discuss a simple way 
to characterize the properties of the collective oscillations 
of a general multi-component mixture. Restricting to 
temperatures larger than the superfluid critical temper- 
ature (see discussion in Section II), we will describe the 
Fermi gas mixture using the Boltzmann equation, with no 
intra-species interaction. Introducing a Gaussian ansatz 
for the distribution function, the equations for the small 
oscillations are determined: a convenient way to clas- 
sify the different regimes in the parameter space is based 
on the study of the eigenvalues of the small oscillations 
matrix. Once that the eigenvalues are found, the exper- 
imental quantities (like the center of mass positions for 
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dipole oscillations) can be determined. The real part of 
the eigenvalues determines the damping of the quantities, 
while the imaginary part is related to the oscillations: we 
show that the number of peaks of the analytical continu- 
ation of the Laplace transform of the physical quantities 
gives the number of real solutions (at variance, the peaks 
of the Fourier transform do not give the number of com- 
plex solutions). One can then classify the solutions of 
the linearized Boltzmann equation through the number 
of real eigenvalues of the small oscillations matrix. A real 
root corresponds to an overdamped eigenmode, whereas 
complex ones correspond to oscillatory modes and purely 
imaginary to undamped oscillatory modes. The small 
oscillation matrix is real, therefore complex eigenvalues 
appear in pairs of complex conjugates and the number 
of oscillation frequencies is simply the number of pairs 
of complex roots. Counting the overdamped eigenmodes 
(i.e. the real roots), one has a measure of the degree of 
collectivity of the system: the larger this number, the 
stronger the effect of interactions and trap confinement 
in reducing the degrees of freedom of the mixture. We 
will show that this approach is complementary to the 
study of the oscillatory frequencies, and convenient to 
classify the properties of the small oscillations solutions. 
We will focus on the case of dipole oscillations (i.e, one- 
dimensional out-of-phase motion of the centers of mass of 
the various components), but the proposed classification 
works for general collective modes of multi-component 
fermionic mixtures. 



The plan of the paper is the following: In the next 
Section we review the formalism of the Boltzmann equa- 
tion and the method of averages, applying it to a multi- 
component fermionic mixture and deriving the equations 
of motion for the centers of mass of the components. 
In Section III we study the dipole oscillations of a two- 
component mixture in harmonic potentials having dif- 
ferent trapping frequencies: indeed, when the trapping 
frequencies are equal, two of the eigenvalues of the 4x4 
small oscillation matrix are purely imaginary due to the 
Kohn theorem [l|, [H, 13 , while for different trapping fre- 
quencies one can have either or 2 or 4 real solutions, 
similarly to what happens for a three-component mixture 
with equal frequencies. In Section III we also recall that 
a convenient way to count the number of real eigenval- 
ues is based on the Sturm theorem and we discuss how 
to characterize the different solutions starting from the 
experimental quantities. In Section IV we study dipole 
oscillations in three-component mixtures, presenting re- 
sults for a mixture of three hyperfine levels of e Li and for 
a mixture of two hyperfine levels of 6 Li and one of 40 K , 
like in the recently reported experimental results [25M27j j . 
In Sec. V we draw our conclusions and discuss possible 
further developments, whereas in the Appendices addi- 
tional material is presented. 



II. BOLTZMANN EQUATION FOR A 
MULTI-COMPONENT MIXTURE 

In this Section we discuss the main properties of 
the Boltzmann equation for a trapped multi-component 
fermionic mixture. Consider a trapped balanced two- 
component Fermi gas; when the temperature T is sig- 
nificantly larger than the Fermi temperature Tp, the 
gas behaves classically and is in the collisionless regime. 
Decreasing T, the frequency of collisions increases lead- 
ing to a hydrodynamical behaviour. Further decreas- 
ing the temperature T (but with T larger than the su- 
perfluid critical temperature), degeneracy effects become 
strong and one would expect the mixture to be colli- 
sionless again. For weak (inter-species) interactions, this 
sequence collisionless-hydrodynamical-collisionless actu- 
ally occurs as T/Tp decreases from values larger than 
one towards smaller ones. However, if the interac- 
tions are strong (i.e., close to the unitary limit) the se- 
quence found in the experiments [22|, [H[ is collisionless- 
hydrodynamical as it would be in a classical gas (i.e, 
without degeneracy effects): the gas remains hydrody- 
namical. This finding can be explained by taking into ac- 
count in-medium effects: the in-medium enhancement of 
the cross-section compensates Pauli blocking at low tem- 
perature leading to a hydrodynamical behaviour [H, [42j . 
In other words, classical statistics seems to work well 
when interactions are strong because Pauli blocking and 
in-medium effects almost cancel each other. In summary, 
depending on the interaction strength, either degeneracy 
dominates or is compensated by medium effects and the 
low temperature behaviour is different in the two cases; 
however, whatever the interaction strength is, for tem- 
peratures above ~ 0.5Tp the Boltzmann equation frame- 
work with no intra-species interactions and classical col- 
lision term with no medium-effects for the inter-species 
interactions is found to work quite well, not only quali- 
tatively, in reproducing the collective modes (23l. |42| . 

We will take advantage of this finding, and describe 
the Fermi nature of the particles simply by forbidding 
intra-species interactions. In multi-component imbal- 
anced mixtures each component has its own Fermi tem- 
perature Tp^ a , therefore we expect our scheme to work 
well above ~ 0.5 max{Tp jQ ,}. 

If we restrict ourselves to the normal phase (that is, for 
temperatures above the critical temperature for superflu- 
idity), each component of the mixture can be described 
by a semiclassical distribution function f a = / Q (r,p,i) 
obeying the Boltzmann equation [451 ] - Here and in the 
following we denote the N components of the multi- 
component mixture by a = 1,---,N, the a-th com- 
ponent having N a atoms of mass m a . The function 
f a (r,p,t) gives the probability of finding a particle of 
species a at time t in the d 3 rd 3 p volume centered in (r, p) 
in phase space. The normalization condition of f a is 

J dTf a {r, p, t) = N a , with dY = d 3 rd 3 p , (1) 
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and the density is 

p a {r,t) = J d 3 pf a (r,p,t) 



(2) 



The Boltzmann equation for our multi-component 
mixture reads 



N 



df a . df a df a _ x v 

dt + a ' 9r + Pa ' 9p ~ ^ 

7=1 



where 



3Va 

9r 



P 



(3) 



(4) 



The right-hand side of Boltzmann equation (J3J) is a sum 
of collision integrals and the ' on the sum means that we 
sum on all 7 ^ a. V a is the single particle potential felt 
by the a species. In general, different components feel 
different anisotropic trapping potentials: 



V Q (r)^(^ 2 



o 2 v 2 



(we will not consider in the following the case in which 
the harmonic potentials V a have different centers j46|V. 
The collision integrals for classical statistics read [45| 



d 3 Pl dn^L 

ail 



P 



Pi 



(/ Q / 7 i-/4/;i) (6) 



and represent the variation of f a due to collisions with 



particles of type 7. In the latter 



is the differential 



cross section of an a atom and a 7 atom in the center 
of mass of the collision and ft is the angle between the 
relative outgoing and relative ingoing momenta of the 
two colliding particles (in our model there is no intra- 
species interaction: da aa = 0). The atoms a and 7 
have respectively momenta p and p x before the scat- 
tering and p' and p^ after: momentum and kinetic en- 
ergy are conserved in the collision [p + Pi = p' + p[ 
and p 2 /2m a + •p\j2m 1 — p' 2 /2m a + pf /2m 7 ] and all 
the / are evaluated at the same point and time, but 
different momenta: f a = / a (r,p,t), / 7l = f 1 (r,p 1 ,t), 
f' a = f a (v,p',t) and/^/^r.p'^). 

With the classical collision term given in Eq. © 
the equilibrium distribution function is the Maxwell- 
Boltzmann distribution 



/a(r,p) 



cxp{-/3[p7(2m a ) + V a (r) - n a ]} 
(2tt/i) 3 



(7) 



where (3 = l/ksT and the chemical potential is fixed by 
the normalization condition Eq. ([1]). For the harmonic 
trapping potential of Eq. (|5|) it is 



/'n 



where 



UJ a = (u> ax U! ay LJ az ) 1 / 3 



(8) 



(9) 



is the trap average frequency. 

The effect of Fermi or Bosc statistics (blocking or anti- 
blocking of the final state of a collision) can be incor- 
porated in the collision integral by adding appropriate 
factors: then, consistently, the equilibrium distribution 
function to which the system is driven by collisions is a 
Fermi or Bose one plj . 

To study collective oscillations without directly solving 
the Boltzmann equation, it is useful to derive equations 
for integrated average quantities: defining the average of 
a generic quantity x — x( r i P) m the component a as 

(x) a = ^J dTf a (r, p, t) X (r, p) , (10) 

from the Boltzmann equation (J3|) one finds 



d (x) a I P _ dx 
dt \ m„ dv 



dr dp 



N 



dV a dx\ __x^i /T \ 



7 = 1 



(5) where we have defined the collisional average 



{tcryX};, = jj- j dTI alX {r,p) 



(11) 



(12) 



The collisional averages of a quantity x satisfy the fol- 
lowing constraint in each pair of components 
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N a {I alX ) a + N 7 (I iaX ) 7 = § / dTd 3 Pl dQ^ 



da n 



P Pi 



where Ax = Xa + X 7 i - x' a ~ X 7 i : summing the latter 
over a, 7 one finds the following global constraint 



a, 7 a, 7 



~dtt 



P Pi 



(/a/71 - f'M A x 



(13) 



If X is a collisional invariant, i.e., if Ax = 451 ]. then 



£' 

a, 7 



N a {I ai x) a = 



(14) 



Eq. (fTTj) gives rise to a set of equations involving differ- 
ent observables: the goal is to have a closed set of equa- 
tions to solve. This approach has been used to determine 
in a classical gas the shift of the collective frequencies due 
to the interparticle collisions in a Bose-Einstein conden- 
sate above the condensation temperature [47]-[5(| and it 
has been applied to the study of the scissor mode in two- 
component Fermi mixtures [391 ] as well as to Fermi-Bosc 
mixtures [5l], [52[ . A convenient way to evaluate the col- 
lisional averages to close the set of equations and study 
the small oscillations for the linearized Boltzmann equa- 
tion is to make a Gaussian ansatz for the time-dependent 
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distribution function /: in [47J the frequency and the 
damping of quadrupole oscillation in a classical gas ob- 
tained using the Gaussian ansatz was compared with the 
numerical solution of the Boltzmann equation showing a 
good agreement. 

For dipole oscillations along the direction x, choosing 
X = x,p x from Eq. (| 1 1 1) one gets the following 2N equa- 
tions 



' d{x) a {p x ) a 



N 



dt m 

d {Px) a 



~ ^ a ~l X ) a 



7=1 N (15) 

a<*>lx ( X )a = - CW^a > 
7=1 



dt 



with a = 1, • • • , N . To close the set of Eqs. (jT5|). adapting 
the method of [47| we make the following Gaussian ansatz 
for the distribution function: 



f a (r,p,t) 



exp 



(P-Po 

2m„ 



V Q (r - r Q ) - fi a 



(27r/l) 3 

(16) 

where p Q = p a (t) and r Q = r a (t) depend on time. 
Eq. (fT5|) is a local equilibrium ansatz: it general- 
izes (JTJ) by giving it a time-dependent average coordi- 
nate [{r) a = r a (t)} and average momentum [(p) Q = 
p Q (t)]; moreover, it preserves the normalization condi- 
tion (J dTf a {r,p,t) =N a ), 

Using the Gaussian ansatz (1161) and linearizing the Eqs. 
(|15|) (i.e., retaining only terms that are at most linear 
in (x) a and (p x ) a ), one finds a closed set of equations. 
Indeed the collisional averages in the right-hand side of 
Eqs. ([15]) can be evaluated giving 



(17a) 



(I ajPx ) a = ^NJ^k-^-) (17b) 

where /i a7 is the two body reduced mass (/i Q7 = 
with m aj = m a + to 7 ): the explicit compu- 
tation is reviewed in Appendix A. In Eq. (|17bl) r Q7 is a 
parameter related to the collisions between atoms a and 
7 and will turn out to be proportional to the relaxation 
time of the dipole mode: r Q7 is defined as 



1 



' a 7 



4/r/r ( oj a uj. 



3?r 



J a7 



7/ ' 



(18) 



where y a7 = „ h ^ .± , with a a7 the scattering length be- 
tween the species a and 7. The function / is defined in 



Eq. (|A4[) and at the unitary limit / = 1; u> a is defined in 
Eq. @ and we also used the notation 



n 

i=x,y,z 



1/6 



(19) 



Notice that repeating the derivation including Fermi 
statistics one would still get the equations Eqs. (JTTJ) , the 



only modification being in the temperature dependence 
of the parameter r Q7 . 

Equations[l5l together with Eqs. (fT7| . form a closed set 
of differential equations describing the small oscillations 
of the centers of mass of a multi-component fermionic 
mixture at T > 0.5 max{JV |Q }. 



A. Two components with equal trapping 
frequencies 

As a simple application of Eqs. flT5"j) we retrieve in this 
Section the well-known case of a two-component mix- 
ture having equal trapping frequencies along the dipole 
motion direction: u)\ x = U2 X = w x . One easily checks 
that the coordinate and velocity of the center of mass 
of the mixture decouple from the relative position and 
velocity and oscillate with frequency ui x , according to 
Kohn theorem. By defining the (average) relative posi- 
tion and velocity respectively as x ic \ : 
Vtei = (Px)i /mi - (p x ) 2 /m 2 one gets 



(a;), — (x) 2 and 



dx T 



v Ic i = 



dt 

aUrol . 2 



-T^rel 

r' 



where the inverse relaxation time is — M1+M2 J_ anc j 

t' mi2 T12 

M a = N a m a is the total mass of the component a. It 
follows that if 



> 2uj x 



(20) 



then the relative motion is overdamped: as long as the 
interactions are larger than a critical value the transport 
of momentum is so effective that the two components 
cannot move independently. 



III. TWO COMPONENTS WITH DIFFERENT 
TRAPPING FREQUENCIES 

In this Section we consider a mixture of two compo- 
nents in potentials having different trapping frequencies 
(u)ix,wiy, ujiz) and (u>2x, <^2y, U2z)- This case deserves a 
separate discussion since the properties of the solutions 
of the small oscillation equations are similar to those of 
a three-component mixture with equal trapping frequen- 
cies: when the trapping frequencies are equal along the 
direction of the dipole motion {lo\ x = ll>2x), two of the 
eigenvalues of the 4x4 small oscillation matrix are purely 
imaginary due to the Kohn theorem [l|, @, E3| , while for 
different trapping frequencies (uji x ^ uj2x) one can have 
either or 2 or 4 real solutions, as it happens for a three- 
component mixture. In this Section we also apply the 
Sturm theorem to count the number of real eigenvalues of 
the small oscillation matrix and we discuss how to char- 
acterize the solutions according to this number starting 
from the experimental quantities. We finally discuss the 
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FIG. 1: Plot of the regions having (dashed) and 2 (blank) 
real roots in the plane a,b (a = — — , b = — 2 — - — ) for 

2 

the case c = ^p- = 4. Inset: plot of regions having (dashed), 
2 (blank) and 4 (dotted) real roots for the case c = 0.001. 




CO s 

FIG. 2: Top: Plot of the center of mass Xi (solid line) and 
X2 (dashed line) as a function of time for a = 1, b = 1 (circle 
in Fig[T|. Bottom left: modulus of the Fourier transform of 
Xi,2, given by (127(1 . Bottom right: modulus of the Laplace 
transform of defined in (|29l) . Initial conditions: Xi(0) = 
1, X 2 (0) = 0, Vi(0) = F 2 (0) = 0; c = 4 as in Figffl 



properties of the Fourier and Laplace transforms of the 
center of mass positions giving informations respectively 
on the imaginary and real part of the eigenvalues. 

We define the dimensionless quantities X a (t) = 
(x) a (t)/L and V Q (i) = (p x ) a (t)/(m a uj lx L) (a = 1,2 and 
L is a length scale of the system which can be set equal 
to an harmonic oscillator length); time will be scaled in 
units of l/coi x . The equations of motion for the dimen- 
sionless centers of mass positions X a and velocities V a 
are then 



X l = V l 

x 2 = v 2 
Vi = -x x - 

V 2 = -cX 2 



b(Vi - 
ha(Vi 



V 2 ) 
-V 2 ) 



where a 



Ah 



b = 



— M 2 



(21) 



'~ : the parame- 



nia TWi x ' mi2 TOl, ' 

ter r is defined in Eq. (ITS)) and M a = N a m a is the total 
mass of the component a. 

It is convenient to introduce a matrix notation for 
the equations of motion: defining the vector Y = 
(Xi, X 2 , Vi,V 2 ) T we can write Eqs. as 



Y = G 2 Y , 
where the 4x4 matrix G 2 is given by 



G 2 = 



( 



-1 

V -c 



1 



-b 
a 



(22) 



(23) 



, A4 of the 



tain the time evolution of the centers of mass. If B 
is the matrix diagonalizing G 2 such that B~ 1 G 2 B = 
diag (Ai, • • • , A4), one has 



Y*(t) 



fc=i 



(24) 



where the coefficients Cjk depend upon the initial condi- 
tions through the relation Cjk = Bjk Yli=i ^/^(O). 

The properties of the solutions of Eqs. (|22|) are com- 
pletely determined by the properties of the eigenvalues 
Afc of the matrix (|23|) , which are given by the roots of the 
equation 

P(A) = A 4 + (a + 6)A 3 + (l + c)A 2 + (a+6c)A + c= , (25) 

where P(A) is the characteristic polynomial of the matrix 
G 2 . Since P(A) is a polynomial with real coefficients, if 
it has a complex root it has also as a root its complex 
conjugate: this means that Eq. ([23)1 has either or 2 or 
4 real solutions. It is easy to verify that the solutions of 
the equation -P(A) = satisfy the property 



\i R) < 



(26) 



Of course, once that the eigenvalues Ai, 

matrix (I23|) have been determined, it is possible to ob- respectively by A^' 1 and A 



which just states the stability of solutions (see Appendix 
B): the real parts of the eigenvalues A& correspond to the 
damping of the normal modes and the imaginary parts to 
their frequencies. In Eq. (|2"rj|) . as well as in the following, 
we denote the real part and the imaginary part of A& 
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One has that = if and only if c = 1 (i.e., u)\ x — 
id 2x ): indeed if c= 1, P(A) = (A 2 +'l)[A 2 + (a + b)X + 1)] 
and two solutions are just ±i. These two roots corre- 
spond to the oscillations of the center of mass of the 
mixture with frequency equal to the trap one and with- 
out damping, as discussed in Section III Al Then for equal 
trapping frequencies Eq. (|25j) can have only or 2 real 
solutions, and the critical value for passing from one re- 
gion to the other is just given by Eq. (|20|) . 

To compute the number of real solutions we could look 
at the general solution of the fourth degree equation (|25|) ; 
however, in order to generalize this approach to three (or 
more) components it is more convenient and straightfor- 
ward to use the Sturm algorithm which we briefly 
recall in Appendix B. In Fig 0] we plot the regions with 0, 
2 and 4 real roots in the plane a, b for two values of c; the 
inset shows that in some narrow region of the parameter 
space (c = 0.001 in the inset) regions with 4 real roots 
may occur. 

If we adopt the complementary point of view of count- 
ing the number of independent oscillations (the number 
of pairs of complex solutions) we see that: only for equal 
trapping frequencies (c = 1) we have an undamped oscil- 
latory mode (Kohn mode); on the contrary, for unequal 
frequencies (c ^ 1), all the modes are either damped 
or overdamped. The number of the damped oscillatory 
modes is 2, 1 or 0, corresponding then to 0, 2 or even 4 
real roots of Eq. (|25|) . 

Fourier and Laplace transforms. We now discuss how 
to determine the regions with a different number of real 
roots from experimentally measured quantities. When 
the centers of mass of the two components have been 
determined or measured, informations on the real and 
imaginary parts of the eigenvalues Xk can be obtained 
respectively from their Fourier and Laplace transforms: 



F a (w) ee / e~ iut X a (t)dt , 



L a (s) 



e- st X a (t)dt 



(27) 



(28) 



Each pair of complex roots, say Ai = A^ + iX^ 

and A2 = A^ — iXi , having an imaginary part HX^\ 
should correspond to a maximum of the modulus of the 
Fourier transform at id = X^ : then, one could expect 
that if the Fourier transform has (for positive values of 
u>) two distinct peaks there should be 4 complex eigenval- 
ues (i.e., real solutions). However, Fourier transforms 
of the centers of mass can be analytically computed for 
arbitrary initial conditions, showing that two peaks (at 
two positive values of id) do not always occur when Eq. 
(j2"5| has 4 complex roots: e.g., if the positions at which 
the two peaks should be are too close, the peaks merge 
into one. Therefore one cannot easily infer the number 
of complex solutions from the peaks of the Fourier trans- 
forms. 



On the other side, to each real solution corresponds a 
pole in the Laplace transform when s is extended to neg- 
ative values: indeed, setting X a (t) = J2t=i Cke Xkt with 
the coefficients Ck depending upon the initial conditions, 
one has for s > the analytical expression 



- (an.) 



4 

£ 

k=l 



Ck 



Xi 



(29) 



Eq. ([2"§]) can be defined also for negative values s and 
it has a pole when s is equal to a real (and negative) 
eigenvalue. 

One can experimentally measure X a (t): of course, to 
determine the eigenvalues Xk one can directly fit the 
data with the expression X a (t) = J2k=iCke Xkt . In 
particular, since the number of real eigenvalues is 0, 2 
or 4, one can use as fitting functions the sum of 4 ex- 
ponentials (corresponding to 4 real eigenvalues), or the 
sum of 2 exponentials and an exponential modulated 
by a sinusoidal (corresponding to 2 real and 2 com- 
plex conjugate eigenvalues), or the sum of 2 exponen- 
tials modulated by sinusoidals (corresponding to 4 com- 
plex eigenvalues, i.e. 2 pairs of complex conjugate eigen- 
values): the best fit among them would determine the 
number of real eigenvalues. Another (possibly comple- 
mentary) method to determine the eigenvalues A^ is to 
extend the definition (l28l) for L a (s) to negative values 
of s. One has that the integral in (|28[) is defined for 

s > s m i n = — min{ A^ |}. Therefore the Laplace 
transform of the center of mass position has a diver- 



gence at s 



if lim. 



i)L a (a) = 0, 



then s m in is the real part of a pair of complex conjugate 
roots (say, Ai and A2) and one has to fit from data the 
corresponding imaginary part A^ (with A2 = — X± ) 
and coefficients C^a! if not, then s m j n is a real root (say, 
Ai) and lim s _,. Smin (s — s m i n )L a {s) = C\. In the latter 
case, one can define a quantity X a (t) = X a (t) — C\e Xlt : 
again, from the divergence of the Laplace transform of 
X a (t) for negative values of s one can determine the next 
eigenvalue [similarly, for a complex root one has to define 
X a (t) ee X a (t) - Cie Alt - C 2 e A2t ]. 

The resulting behaviour of the centers of mass posi- 
tions and of the modulus of the Fourier transform (|27|) 
and the Laplace transform (f2"9"| is plotted in FigJ5] and 
Figj3]for two sets of parameters, one belonging to the re- 
gion of zero and one in the region of two real eigenvalues: 
the Laplace transforms show respectively zero and two 
sharp peaks. 



IV. THREE COMPONENTS 

In this Section we discuss the properties of the modes of 
three-component mixtures: we consider three-component 



fermionic mixtures of Li, as in the setup reported in 25] , 



and mixtures with two Li and one K species [2 
We focus for sake of simplicity on isotropic potentials 



7 



V a with the same trapping frequency uj. Rescaling the 
variables as in Section III, measuring time in units of 1/w 
and defining the vector Y = (X 1: X 2 , X 3l Vl, V%, Vz) T we 
have Y = G3Y, where the 6x6 matrix G3 is: 

/ \ 








i 






— Ti2 — IY3 


Tl2 


Tl3 


-1 


r 2 i 


— T21 — T23 


T23 




r 3 i 


1^32 


-r 3X - 



(30) 



where and 1 are the 3x3 zero and identity matrices 
and T a7 ee ^-^t^(1-<5q 7 ) (a, 7 = 1,2,3) with S aj the 
Kronecker delta. The six non-zero parameters r a7 are 
not independent: they satisfy the relation T^I^rsi = 
ri3r3 2 r2i- The characteristic polynomial is: 



P(A) = (A 2 + 1) 



A 4 + A 3 ^r Q7 + A 2 (2+ r « 



a, 7 



E 

a,7,e^a 



r r 

x Q7 1 67 



A$> 



(31) 



The discussion proceeds now as in Section III: one has 



• Aj^ < (Stability of solutions) 



A 



(R) 



& A 



±1 



the number of real solutions can be either 0, 2 or 4 and 
the number of real roots can be inferred from the Laplace 
transforms of experimental measured quantities. At the 
transition line between one region with or 2 real eigen- 
values and a region having more real roots, S functions 
appear and the corresponding strengths of these 5 func- 
tions become non-vanishing. The eigenvalues A = ±i 
correspond to the center of mass oscillation with the fre- 
quency of the trap (Kohn theorem). Again, with differ- 
ent trapping frequencies, regions with 6 real eigenvalues 
(that is, with no oscillatory mode at all) may occur. 



A. Three-component mixture of the same species 

We discuss in this Section the behaviour of the eigen- 
values of a mixture in which the components belong to 
the same species, as in the three-component 6 Li mixture 
described in [25j . 

For an unpolarized mixture in which each species con- 
tains an equal number of particles Nn (SN^i is the total 
number of particles in the mixture), there are three in- 



dependent parameters T Q7 
matrix G3 reads 



(since r Q7 = T ia ) and the 



G, 



-(6+6) 

6 
6 



6 

-(6+6) 
6 



6 
6 

-(6 + 6)/ 




10 



10 



, 



CO 



(32) 



FIG. 3: Top: Plot of the center of mass X\ (solid line) and X2 
(dashed line) as a function of time for a = 3, b — 3 (triangle 
in FigfTJ). Bottom left: modulus of the Fourier transform 
of Xi_2- Bottom right: modulus of the Laplace transform of 
Xi, 2 defined in (J29J). Initial conditions: Xi(0) = 1, X 2 (0) = 0, 
Vi(0) = V 2 (0) = 0; c = 4 as in Figffl 



where £1 = Fi2, 6 = Ti3 an d 6 = -^23- The character- 
istic equation and its solutions are studied in Appendix 
B. 

In Fig |31 we plot the results of our analysis with the 
Sturm theorem in the plane 6"6 f° r two values of £3. If 
£3 < 1 , then one has three different regions corresponding 
to 0, 2 and 4 real solutions. When £3 > 1 the region 
with real solutions disappears and the motion is ruled 
either by one or two oscillating frequencies, but not three. 
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FIG. 4: Plot in the £1-^2 plane of regions having (dashed), 
2 (blank) and 4 (dotted) real roots for a three-component 
mixture having £3 = 0.2 (left) and £3 = 1.5 (right) - the value 
of £3 at which the region with real roots disappears is £3 = 1. 
The points denoted in the left plot are at £1 = 0.5, £2 = 0.5 
(circle), £1 = 1.5, £2 = 1.5 (triangle), £1 = 1.8, £2 = 1.7 
(square) . 



In Figj5] the center of mass positions for a simple initial 
condition and the modulus of their Fourier and Laplace 
(|29f transforms are plotted for three sets of parameters, 
belonging to regions with respectively 0, 2 and 4 real 
eigenvalues showing that the Laplace transforms display 
respectively zero, two and four sharp peaks. 

When the mixture is polarized, one has to use Eq. ([3TH) : 
varying the number of the atoms in the components, it 
is possible to explore the different regions through the 
polarization. An example is shown in FiglHl where the 
number of atoms of the component 3 is varied and the 
components 1 and 2 are at the unitary limit. As dis- 
cussed in [2^| , the mixture is collisionally stable for small 
values of the other scattering lengths. Fig(B] shows that 
by varying the trapping frequency one can explore the 
three regions with 0, 2 and 4 real eigenvalues. Notice 
that the stronger the confinement due to the trap, the 
fewer independent oscillatory modes are possible. 



B. Different species 

In this Section we present the results for a three- 
component mixture of two species of 6 Li (denoted by 
Lii and H2) and one species of 40 K recently realized 
experimentally: a small sample of 40 K is immersed in a 
two-component 6 Li mixture and found to be collisionally 
stable with very low loss rates for large negative scatter- 
ing lengths between the two 6 Li (for | an-u |> 1400a , 
where ao is the Bohr radius) |27j ■ The matrix G3 of Eq. 



FIG. 5: Plot of the center of mass Xifi,z(t) (top row), modulus 
of the Fourier transform ^1,2,3(0?) (central row) and modulus 
of the Laplace transform £1,2,3(5) defined in (|29[) (bottom 
row). The left, central and right columns correspond respec- 
tively to the point £1 = 0.5, £2 = 0.5 (circle in Fig(4]) - real 
roots), £i = 1.5, £2 = 1.5 (triangle - 2 real roots), £1 = 1.8, 
£2 = 1.7 (square - 4 real roots). Solid, dashed and dotted 
lines are for the corresponding quantities of the species 1, 2 
and 3. Initial conditions: Xi(0) = 1, X 2 (0) = 0.5, X 3 (0) = 0, 
Vi(0) = Va(0) = 14(0) = 0. 



(|3TJ)) simplifies to 



q 3 - 

£1 -(£i+6™0 &ra 

\ 6 6 -26/ 

(33) 

where we defined £1 = Ti2, £2 = r = mx/wiLi and 
a = Nk/Nh. To simplify formulas, we assumed the 
same number of atoms in the two Li species and equal 
isotropic trapping frequencies, but similar (even if more 
involved) expressions can be found for the realistic case 
of anisotropic cigar-shaped potentials and with the trap- 
ping frequencies of Li\-Li2 different from those of the K 
component: they will be reported elsewhere. 

In Fig[7]we plot the number of real solutions taking a 
mixture with the same number of particles for the two 
Lithium components and a smaller number of particles of 
the K component. We vary the scattering length between 
Li\ and Lii and the trapping frequency u>. 

The number of oscillatory modes is reduced, again, by 
increasing the interaction strength or making the trap 
more conhning. 
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FIG. 6: Plot of regions having (dashed), 2 (blank) and 4 
(dotted) real roots for a mixture of three Li species when the 
number N3 and the trapping frequency u = 2nv axe varied. 
The components 1 and 2 are taken at unitary limit with a 
fixed number of atoms Ni = N% — 5 ■ 10 4 . Temperature is 
fixed (T = 500nK) and the two remaining scattering lengths 
are chosen to be 013 = 023 = 400ao (similar plots are obtained 
when O13 and 023 are different, but small). 




FIG. 7: Plot of regions having (dashed) and 2 (blank) real 
roots for a mixture of Li\-Hi-K with equal trapping fre- 
quencies uj — 2ttv and equal number of Li atoms Nn x = 
Nhi 2 = Nli = 5 ■ 10 4 . The number of K atoms has 
been chosen Nk = 0.25NLi- Parameters: T = 200nK, 
a,L%x—K = a,Li 2 -K = 63ao; the Li\ — Lii scattering length 
ai2 is measured in units of lOOao. For v < 55Hz there are 
real roots whatever 012 is. 



V. CONCLUSIONS 

In this paper we used the Boltzmann equation to study 
dipole oscillations in a general fermionic mixture: to dis- 
cuss and classify the different solutions in the parameter 
space through it, we used the number of real eigenval- 
ues of the small oscillations matrix, discussing how this 
number is related to the Laplace transform of the exper- 
imental quantities. The small oscillation matrix for a N- 
component mixture has 27V eigenvalues and the complex 
ones appear in pairs of complex conjugates: the number 
of independent collective mode frequencies is therefore 
given by (2N — number of real roots)/2. Counting the 
real roots (overdamped modes) and the complex ones 
(oscillatory ones) are complementary descriptions. The 
reason why we chose the first one is to take advantage of 
Sturm algorithm and Laplace transform properties. 

After deriving the equations of motion for the centers 
of mass of the components and discussing as a simple ex- 
ample the case of two components in harmonic potentials 
having different trapping frequencies, we studied dipole 
oscillations in three-component mixtures. 

A Af-component mixture where each component expe- 
riences a different trapping has a maximum of N inde- 
pendent oscillation frequencies. We found, as expected, 
that when interactions are favoured (by increasing the 
strength of interactions between components, or shrink- 
ing the traps) the components tend to move together and 
the number of independent frequencies reduces. If all 



the trap frequencies are equal along one direction there 
is one undamped oscillation with frequency equal to the 
trap one (Kohn theorem). Beside that, there can be up 
to N — 1 damped oscillatory modes. On the contrary, if 
the trap frequencies are different, there is no purely oscil- 
latory mode. The number of damped oscillatory modes 
can be at most A^ and can be reduced down to 0: it is 
possible that there is no collective dipole oscillatory mode 
at all in the mixture. 

In explicit computations we used the Boltzmann equa- 
tion with classical collision terms and without intra- 
species interactions: this gives, for strong inter-species 
interactions, results very close to the findings obtained 
including Fermi statistics and in-medium effects thanks 
to a compensation of effects [23|, |42| ■ Using the classical 
Boltzmann equation without intra-species interactions is 
expected to be a good approximation (also far from the 
unitary limit) at temperatures above ~ 0.5 max{TF jCt }, 
where Tf, u is the Fermi temperature of the component 
a: for two components, e.g., it well reproduces the ex- 
perimental results for the scissor mode for these tempera- 
tures [H, |42j ■ We also observe that the explicit inclusion 
of the Fermi statistics modifies the dependence of the re- 
laxation times upon the system parameters, but not the 
dynamical equations for the centers of mass and their 
momenta and then the subsequent analysis based on the 
Sturm theorem would proceed in the same way. 

The focus of this paper has been on the study of dipole 
oscillations, however we point out that a similar study 
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can be done for different (and more complicated) collec- 
tive modes: one has to write the small oscillations equa- 
tions, and determine the number of real roots of the small 
oscillations matrix. 

To conclude, we observe that multi-component mix- 
tures are promising to study mechanisms of color super- 
conductivity [33j due to the possibility of having many 
components with tunable interactions: from this point 
of view, to study the collective oscillations (and follow 
the modifications and the fate of the regions in the pa- 
rameter space discussed in this paper) when the temper- 
ature is lowered until the superfluid critical temperature 
can provide useful informations on possible intermediate, 
new strongly correlated phases. 
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FIG. 8: Plot of the function f(y) denned in Eq. (|A4[) as a 
function of the dimensionless parameter y = ti 2 /3/2fj,ai2- 



Appendix A: Computation of r 

In this appendix we give some details on the computa- 
tion of the parameter r related to the collisions between 
atoms of species 1 and atoms of species 2. As usual, 
given two particles with masses 7771,777,2, positions ri,r2 
and momenta Pi,P2 the global and relative coordinate 
and momenta are defined as: P = p 1 +p 2 , R = (toiTi + 
m 2 r 2 )/mi2, r rcl = n - r 2 , p rcl = (m 2 P! - mip 2 )/m 12 . 
Then pj/2m 1 + p 2 /2m 2 = P 2 /2m 12 + p 2 cl /2fi and p. 
/ir rc i, where /_t 



rel 



is the reduced mass and m\ 2 



mi+m,2 

mi +m2. For simplicity we consider equal isotropic trap- 
ping frequencies: wj.^ = uj 2 .i = lu (i = x,y,z); a\ 2 is the 
scattering length and ^ = p s /tp+i/a 1 * nc differential 
cross section. The collisional average (Ii 2 p x )i is defined 
as 



{I\2Px) 1 = jjr J ' dTI l2 p x 



(Al) 



Plugging the Gaussian ansatz for the distribution func- 
tion (fT6|) in Eq. (|Al|) and moving to the center of mass 
and relative coordinates we obtain 



{Ii2Px)x = : — J aPrei-^r 



3tt 



Prel 



-e 2 f* 



r,-' 



(A2) 

where fi is the reduced mass and v vc i the relative velocity 
between the two components. Setting j^P^ ci = t then we 
have: 



IT . m 2 p 2 uj\ T r°° , 

3vr J t + y 



(A3) 



where y 



h 2 /3 



The last integral in the right-hand 
side of (|A3|) can be written in terms of the generalized 
T-function r(0,x) = f™dt^ H as 



Jo t + y 



e- f = l-2/ + yVr(0,y) . (A4) 



The function f(y) is plotted in FigU it is /(0) = 1 (at 
the unitary limit) and / vanishing for large y. One then 
obtains 



^ 2 h 2 us- 

3tt 



f(y) 



(A5) 



The computation for general anisotropic and different 
trapping frequencies can be similarly done and gives the 
result (fl8l). 



Appendix B: Sturm chains and Sturm theorem. 

Here we recall how to apply Sturm theorem to count 
the number of real roots of a polynomial equation P(A) = 
53]. First, one has to define the (finite) sequence of 
polynomials (Sturm chain) as follows: 

Po(X) = P(A) 

P 2 (A) = -rem(P (A),Pi(A)) 
P 3 (A) = -rcm(P 1 (A),P 2 (A)) 

= -rem(P m _i(A),P m (A)) , 

(Bl) 
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where we denote with rem(P m _i(A), P m (X)) the rest of 
the polynomial division of the polynomials P m _i(A) and 
P m (A). Then one has to evaluate the sign of the polyno- 
mials Pj (A) at +00 and —00 and to compute the number 
of changes of sign (which we call A(+oo) and A(— 00)) 
in the sequences {Po(+oo), Pi (+00), P2(+oo), . . . } and 
{Po(— 00), Pl(— 00), P 2 (— 00), . . . }. Sturm theorem states 
that the number of real solutions of the equation P(x) = 
is given by |A(+oo) — A(— 00) |. Notice that the Sturm 
theorem is not directly applicable to cases where there 
are degenerate solutions. In these cases, one has first 
to transform the equation into one that has only simple 
roots, but this can always be done. 

We show how to apply Sturm theorem in the specific 
case of polynomial (|25[) of fourth degree 

P(A) = A 4 + A 3 (a + 6) + A 2 (l + c) + A(a + 6c)+c . (B2) 

In order to evaluate the sign of Sturm polynomials at 
A = ±00 we determine the sign of the coefficient of the 
highest power in the polynomial as follows: 



A 


Po Pi P2 


P 3 Pi 


+00 


+ + (-1)" 


(-1)" (-1) 7 


—00 


+ - (-1)" 


-(-IF (-1) 7 



where a, /3, 7 = {0,1}. Then number of changes of signs 
in the two cases is: 

A(+oo) = |_i[(-l)° + (-l)°(-l)^ + (-l)P(-l)7] 
A(-oo) = ! + l [(-!)" + (-l)°(-l)^ + (-l)P(-l)7] 
and the function counting the number of real solutions is 

|A(+ ( x 5 )-A(-oo)| = |l + (-l)° + (-l)°(-l)' 3 + (-l)^(-l)-'| . (B3) 

The counting function (|B3|) is then 





a=0 


a=l 




7=0 7=1 


7=0 7=1 


/3=0 


4 2 


2 


0=1 


2 


2 



When the real part of the root of Eq. (|B2I) is negative, 
then solutions of equations of motion (|2"Tj) do not grow 
exponentially: writing P(A) = A 4 + AX 3 + BX 2 + CX + D 
this happens if [HI 

A>0 P > C>0 P> > 
ABC > A 2 D + C 2 . 

It is easy to verify that these conditions arc satisfied. 

The characteristic polynomial of the matrix G3 in 
Eq.® is: 

P(A) =(A 2 + 1){A 4 + 2A 3 (6 + 6 + 6) 

+ A 2 [2 + 3(a6+66 + ae 3 )] 

+ 2A(6+6+6) + l}- 
For completeness we write here the explicit solutions 
of the quartic equation P(A)/(A 2 + 1) = 0: 

Xiji = i (-S- v /Qi±^/-4 + Q 2 + 2S v / Qi) 
X111.1v = \ (-S-y/Q^±^-4: + Q 2 -2S^/Q^j , 

where 

5(^,6,6) = 6+6 + 6 
Qi (6 , 6, 6) = 6 2 + 6 2 + 6 2 - 66 - 66 - 66 
Q 2 (6,6,6) = 2(6 +6 + 6) 2 - 3(66 + 66 + 66) ■ 



[1] L. P. Pitaevskii and S. Stringari, Bose-Einstein conden- 
sation (Oxford, Oxford University Press, 2003). 

[2] C. J. Pethick and H. Smith, Bose-Einstein condensation 
in dilute gases, 2nd ed. (Cambridge, Cambridge Univer- 
sity Press, 2008). 

[3] W. Ketterle and M. W. Zwierlein, in Ultracold Fermi 
Gases, Proceedings of the international school held at 
Villa Monastero, Varenna, 20-30 June, 2006 eds. M. In- 
guscio, W. Ketterle, and C. Salomon (Amsterdam, IOS 
Press, 2007). 

[4] S. Giorgini, L. P. Pitaevskii, and S. Stringari, Rev. Mod. 

Phys. 80, 1215 (2008). 
[5] B.DeMarco and D.S. Jin, Science 285, 1703 (1999). 
[6] M. Greiner, C. A. Regal, and D. S. Jin, Nature 426, 537 

(2003). 

[7] S. Jochim, M. Bartenstein, A. Altmeyer, G. Hendl, C. 
Chin, J. H. Denschlag, and R. Grimm, Science 302, 2101 



(2003) . 

[8] M. W. Zwierlein, C. A. Stan, C. H. Schunck, S. M. F. 

Raupach, S. Gupta, Z. Hadzibabic, and W. Ketterle, 

Phys. Rev. Lett. 91, 250401 (2003). 
[9] T. Bourdel, L. Khaykovich, J. Cubizolles, J. Zhang, 

F. Chevy, M. Teichmann, L. Tarruell, S. J. J. M. F. 

Kokkelmans, and C. Salomon, Phys. Rev. Lett 93, 050401 

(2004) . 

[10] G. B. Partridge, K. E. Strecker, R. I. Kamar, M. W. Jack, 

and R. G. Hulet, Phys. Rev. Lett 95, 020404 (2005). 
[11] J. Kinast, A. Turlapov, J. E. Thomas, Q. Chen, J. Stajic, 

and K. Levin, Science 307, 1296 (2005). 
[12] G. Veeravalli, E. Kuhnle, P. Dyke, and C. J. Vale, Phys. 

Rev. Lett. 101, 250403 (2008) 
[13] G. Modugno, F. Ferlaino, R. Heidemann, G. Roati, and 

M. Inguscio, Phys. Rev. A 68, 011601(R) (2003). 
[14] J. K. Chin, D. E. Miller, Y. Liu, C. Stan, W. Setiawan, C. 



12 



Sanner, K. Xu, and W. Ketterle, Nature 443, 961 (2006). 
T. Stoferle, H. Moritz, K. Giinter, M. Kohl, and T. [35 
Esslinger, Phys. Rev. Lett. 96, 030401 (2006). [36 
U. Schneider, L. Hackermuller, S. Will, Th. Best, I. [37 
Bloch, T. A. Costi, R. W. Helmes, D. Rasch, and A. 
Rosch, Science 322, 1520 (2008). [38 
M. W. Zwierlein, A. Schirotzek, C. H. Schunck and W. 
Ketterle, Science 311, 492 (2006) [39 
G. B. Partridge, W. Li, R. I. Kamar, Y. A. Liao, R. G. 
Hulet, Science 311, 503 (2006). [40 
S. D. Gensemer and D. S. Jin, Phys. Rev. Lett. 87, 173201 
(2001). [41 
J. Kinast, S. L. Hemmer, M. E. Gehm, A. Turlapov, and 
J. E. Thomas, Phys. Rev. Lett. 92, 150402 (2004). [42 
M. Bartenstein, A. Altmeyer, S. Riedl, S. Jochim, C. 
Chin, J. H. Denschlag, and R. Grimm, Phys. Rev. Lett. [43 
92, 203201 (2004). [44 
M. J. Wright, S. Riedl, A. Altmeyer, C. Kohstall, E. R. [45 
S. Guajardo, J. H. Denschlag, and R. Grimm, 
Phys. Rev. Lett. 99, 150403 (2007). 

S. Riedl, E. R. SanchezGuajardo, C. Kohstall, A. Alt- [46 
meyer, M. J. Wright, J. H. Denschlag, R. Grimm, G. M. 
Bruun, and H. Smith, Phys. Rev A 78, 053609 (2008). [47 
S. Nascimbene, N. Navon, K. J. Jiang, L. Tarruell, M. 
Teichmann, J. McKeever, F. Chevy, and C. Salomon, [48 
Phys. Rev. Lett. 103, 170402 (2009). 

T. B. Ottenstein, T. Lompe, M. Kohnen, A. N. Wenz, [49 
and S. Jochim, Phys. Rev. Lett. 101, 203202 (2008). 

E. Wille et al, Phys. Rev. Lett. 100, 053201 (2008). [50 

F. M. Spiegelhalder, A. Trenkwalder, D. Naik, G. Hendl, 

F. Schreck, and R. Grimm, Phys. Rev. Lett. 103, 223203 [51 
(2009). 

A. Rapp, G. Zarand, C. Honerkamp, and W. Hofstetter, 
Phys. Rev. Lett. 98, 160405 (2007). [52 
S. Capponi, G. Roux, P. Lecheminant, P. Azaria, E. 
Boulat, and S. R. White, Phys. Rev. A 77, 013624 (2008). 
X. W. Guan, M. T. Batchelor, C. Lee, and H.-Q. Zhou, [53 
Phys. Rev. Lett. 100, 200401 (2008). 

T. N. De Silva, Phys. Rev. A 80, 013620 (2009). [54 
P. F. Bedaquea and J. P. D'Incao, Ann. Phys. 324, 1763 
(2009). 

M. G. Alford, A. Schmitt, K. Rajagopal, and T. Schafer, [55 
Rev. Mod. Phys. 80, 1455 (2008). 

K. Rajagopal and F. Wilczek, At the frontier of par- 
ticle physics, World Scientific, Vol. 1-3, Pages: 2062- 



2151,(2001); | arXiv.hep-ph/001 1 333 , 

L. Vichi and S. Stringari, Phys. Rev. A 60, 4734 (1999). 

L. Vichi, J. Low Temp. Phys. 121, 177 (2000). 

G. M. Bruun and B. R. Mottelson, Phys. Rev. Lett. 87, 

270403 (2001). 

M. Tomoyuki and G. F. Bertsch, Phys. Rev. A 73, 013610 

(2006) . 

G. M. Bruun and H. Smith, Phys. Rev. A 76, 045602 

(2007) . 

A. Lazarides and B. Van Schaeybroeck, Phys. Rev. A 77, 
041602(R) (2008). 

H. P. Dahal, S. Gaudio, J. D. Feldmann, andK. S. Bedell, 
Phys. Rev. A 78, 035601 (2008). 

S. Chiacchiera, T. Lepers, D. Davesne, and M. Urban, 

Phys. Rev. A 79, 033613 (2009). 

G. Catelani and E. A. Yuzbashvan. \arXiv:0909.3851\ 
W. Kohn, Phys. Rev. 123, 1242 (1961). 
G. E. Uhlenbeck, G. E. Ford, and E. W. Montroll, Lec- 
tures in statistical mechanics (Providence, RI, American 
Mathematical Society, 1963). 

J. Williams, R. Walser, J. Cooper, E. A. Cornell, and M. 

Holland, Phys. Rev. A 61, 033612 (2000). 

D. Guery-Odelin, F. Zambelli, J. Dalibard, and S. 

Stringari, Phys. Rev. A 60, 4851 (1999). 

U. Al Khawaja, C. J. Pethick, and H. Smith, J. Low 

Temp. Phys. 118, 127 (2000). 

P. Pedri, D. Guery-Odelin, and S. Stringari, Phys. Rev. 
A 68, 043608 (2003). 

W. J. Mullin and R. J. Ragan, Phys. Rev. A 74, 043607 
(2006). 

G. Ferrari, M. Inguscio, W. Jastrzebski, G. Modugno, 
G. Roati, and A. Simoni, Phys. Rev. Lett. 89, 053202 
(2002). 

F. Ferlaino, R. J. Brecha, P. Hannaford, F. Riboli, G. 
Roati, G. Modugno, and M. Inguscio, J. Opt. B: Quan- 
tum Semiclassical Opt. 5, S3 (2003). 

S. Basu, R. Pollack, and M. Roy, Algorithms in real al- 
gebraic geometry (Berlin, Springer, 2003). 
M. Abramowitz, I. A. Stegun, Handbook of mathematical 
functions with formulas, graphs, and mathematical tables 
(Dover, 1972, repr. ed). 

G. Birkhoff and S. MacLane, A survey of modern algebra 
(Macmillan, New York, 1953). 



